Семинары 1-3. Эксперименты.

0.1 0. Организационные моменты

• Мои контакты
  • Почта – annastavnychuk@gmail.com
  • Telegram – https://t.me/AnnaStavniychuk
• Ресурсы по курсу
  • Страница курса на on.econ – https://on.econ.msu.ru/course/view.php?id=1624
  • Чат курса – https://t.me/Applied_ECM_2021_chat
  • Канал курса –
• Не стесняйтесь задавать свои вопросы в чат курса
  • Ответ на ваш вопрос может быть полезен вашим однокурсникам
  • В чате сразу четыре преподавателя, что в несколько раз повышает вероятность получить ответ на ваш вопрос скорее :)
  • Вы также можете отвечать на вопросы ваших однокурсников или дискутировать с ними в чате

0.2 1. Установка R

Установка R на компьютер состоит из двух частей:

• Установка дистрибутива R – https://cran.r-project.org/
• Установка графического интерфейса. Они бывают разными, вы можете выбрать любой удобный для вас. Мы на занятиях будем использовать RStudio – https://www.rstudio.com/products/rstudio/download/

Более подробную инструкцию для разных ОС можно почитать тут.

0.3 2. Короткое напоминание теории

0.3.1 2.1. Потенциальные исходы

Обозначения:

• \(X_i\) – независимые переменные (covariates)

• \(T_i\) – бинарная переменная воздействия (treatment variable): \[\begin{equation*} T_i = \begin{cases} 1, &\text{воздействие на объект i оказано}\\ 0, &\text{воздействие на объект i не оказано} \end{cases} \end{equation*}\]

• \(Y_{i1}\), \(Y_{i0}\) – потенциальные исходы (potential outcomes)

Наблюдаемые исходы \(Y_i\) отличаются от потенциальных исходов. Потенциальные исходы являются гипотетическими случайными величинами, когда наблюдаемые исходы являются фактическими случайными величинами. Наблюдаемые исходы являются функцией от потенциальных исходов:

\(Y_i = T_i \cdot Y_{i1} + (1-T_i) \cdot Y_{i0}\)

Действительно при \(T=1\) мы получим \(Y_i = Y_{i1}\), а при \(T=0\) получим \(Y_i = Y_{i0}\).

Тогда наш эффект воздействия для конкретного наблюдения равен разнице между двумя состояниями мира для этого наблюдения (потенциальными исходами):

\(\tau_i = Y_{i1} - Y_{i0}\)

Фундаментальная проблема причинного вывода (Fundamental problem of causal inference):

Чтобы оценить эффект воздействия для конкретного индивида, мы должны знать потенциальные исходы сразу для двух его состояний мира, но реально мы наблюдаем только одно из них – либо \(Y_{i1}\), если индивид подвергся воздействию, либо \(Y_{i0}\), если он ему не подвергался. Оценка индивидуального эффекта требует доступа к данным, которых у нас физически не может быть.

0.3.2 2.2. Средние эффекты

Если с распределением индивидуального эффекта воздействия (treatment effect) работать не получается, будем довольствоваться средними величинами. Например, попробуем рассчитать средний эффект воздействия (average treatment effect):

\(ATE = \mathbb{E}[\tau_i] = \mathbb{E}[Y_{i1} - Y_{i0}] = \mathbb{E}[Y_{i1}] - \mathbb{E}[Y_{i0}] \xrightarrow{p} \frac{1}{N_1}\sum \limits_{i=1}^{N_1}Y_{i1} - \frac{1}{N_0}\sum \limits_{i=1}^{N_0}Y_{i0}\)

Следующий шаг, который мы предпримем, попробуем рассчиать величину эффекта только для тритмент группы – средний эффект воздействия на задействованных (average treatment effect for the treatment group):

\(ATT = \mathbb{E}[\tau_i|T_i=1] = \mathbb{E}[Y_{i1} - Y_{i0}|T_i=1] = \mathbb{E}[Y_{i1}|T_i=1] - \mathbb{E}[Y_{i0}|T_i=1]\)

А теперь то же самое, но для контрольной группы – средний эффект воздействия на незадействованных (average treatment on the non-treated):

\(ATnT = \mathbb{E}[\tau_i|T_i=0] = \mathbb{E}[Y_{i1} - Y_{i0}|T_i=0] = \mathbb{E}[Y_{i1}|T_i=0] - \mathbb{E}[Y_{i0}|T_i=0]\)

Обратите внимание, как и в случае с определением эффекта воздействия на индивидуальном уровне, различные модификации средних эффектов воздействия снова требуют от нас знания обоих потенциальных исходов для каждого наблюдения. Таким образом, и средние, и индивидуальный эффект воздействия нельзя напрямую рассчиать, но мы будем пробовать их оценить.

Самая простая идея для оценки ATE, которая всем придет в голову, взять простую разницу в средних: \(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0]\).

Но тут всё не так просто, после небольших преобразований, с которыми можно ознакомиться в части 4.1.3 учебника мы получим следующее:

\(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0] = \underbrace{\mathbb{E}[Y_1] - \mathbb{E}[Y_0]}_{\text{ATE}} + \underbrace{\mathbb{E}[Y_0|T=1] - \mathbb{E}[Y_0|T=0]}_{\text{Selection Bias}} + \underbrace{(1-\pi)(ATT - ATnT)}_{\text{Heterogeneous treatment effect bias}}\)

• ATE – интересующий нас эффект
• Selection Bias – смещение, возникающее из-за того, что контрольная группа и группа воздействия различались, даже если бы на них не было оказано воздействие, то есть имеет место некоторый дисбаланс
• Heterogeneous treatment effect bias – различие в интенсивности эффекта для тритмент и контрольной группы, взвешенное на долю выборки \((1-\pi)\), которая попала в контрольную группу

Далее мы рассмотрим предпосылки, которые позволяют нивелировать влияние этих двух смещений и получить оценку ATE.

0.3.3 2.3. Предпосылки

1. Экзогенность воздействия (Independence assumption)

Экзогенность воздействия (Independence assumption) означает, что распределение объекта в тритмент или контрольную группы осуществляется случайно и независимо от его изначальных характеристик. Данная предпосылка обычно обозначается следующим образом \((T_1, Y_0, X)_i \perp T_i\)

Технически для нас это значит следующее:

• \(\mathbb{E}[Y_0|T=1] - \mathbb{E}[Y_0|T=0] = 0 \Rightarrow \text{Selection Bias}=0\)
• \(\mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_1|T=0] = 0\)
  • \((1-\pi)(ATT - ATnT) = (1-\pi)\left[(\mathbb{E}[Y_1|T=1]-\mathbb{E}[Y_0|T=1])-(\mathbb{E}[Y_1|T=0]-\mathbb{E}[Y_0|T=0])\right] = 0 \Rightarrow \text{Heterogeneous treatment effect bias}=0\)

То есть хорошая рандомизация, а следовательно, и выполнение предпосылок, позволяет нам очистить эффект воздействия от двух типов смещения, в этом случае:

\(ATE = \mathbb{E}[Y_1] - \mathbb{E}[Y_0] = \mathbb{E}[Y_1|T=1] - \mathbb{E}[Y_0|T=0] \xrightarrow{p} \frac{1}{N_1}\sum \limits_{i=1}^{N_1}Y_{i1} - \frac{1}{N_0}\sum \limits_{i=1}^{N_0}Y_{i0}\)

2. Отсутствие “внешних эффектов” воздействия (SUTVA – Stable unit treatment value assumption)

Эта предпосылка подразумевает выполнение двух вещей. Во-первых, воздействие оказывается только на один объект и внешние эффекты у него отсутствуют. Во-вторых, воздействие гомогенно – существует только один тип тритмента.

0.4 3. Работа со случайными числами в R

В рамках курса нам неоднократно придется прибегнуть к генерации случайных чисел из какого-нибудь закона распределения, как правило, из нормального. Для работы со случайными числами из нормального распределения в R используется несколько функций: dnorm() для плотности вероятности, pnorm() для функции распределения, qnorm() для квантилей распределения и rnorm() для генерации случайных чисел. Почитать подробнее про них и посмотреть примеры можно тут или тут.

Если вы не уверены в синтаксисе какой-то функции, вы можете воспользоваться встроенной справкой в RStudio (справочная информация появится во вкладке Help):

?rnorm

Если вам когда-нибудь понадобятся функции для других распределений, вы можете, конечно, их просто загуглить или ввести в окне Help запрос “Distributions”:

?Distributions

Теперь, когда мы точно знаем, как устроены аргументы у функции rnorm(n;mean;sd), давайте попробуем достать 5 чисел из стандартного нормального распределения \(\mathcal{N} \sim \left(0;1\right)\):

rnorm(5,0,1)

## [1]  0.25331851 -0.02854676 -0.04287046  1.36860228 -0.22577099

Если вы повторяли за мной, то у вас, вероятно, получились другие значения. А теперь я тоже попробую ещё раз:

rnorm(5,0,1)

## [1]  1.5164706 -1.5487528  0.5846137  0.1238542  0.2159416

Как и ожидалось, у меня тоже не совпало с предыдущим результатом. Существует по меньшей мере две причины, почему нам хотелось бы получать одинаковые результаты:

• Работа на семинарах, когда нам с вами было бы удобно сверяться
• Воспроизводимые научные исследования, когда любой читатель может реплицировать и проверить на корректность ваш результат

В R существует несколько алгоритмов, позволящих генерировать случайные числа (Random Number Generator (RNG), подробнее – ?RNGkind). По умолчанию используется Mersenne twister (Вихрь Мерсенна), поскольку он работает наиболее быстро. Однако все эти алгоритмы на самом деле детерминированы, поэтому генерируемые ими числа корректнее называть “псевдослучайными”.

Генератор псевдослучайных чисел начинает свою работу с определенной точки в пространстве возможных чисел. Эту точку приянто называть начальное число или seed. Начальное число – это число (или вектор), используемое для инициализации генератора псевдослучайных чисел. Если вы хотите получить одинаковые результаты с помощью генератора случайных чисел, важно установить начальное число. Для этого в R используется функция set.seed(). По умолчанию начальное число не установлено, если ничего не указать, то создается новое из текущего времени и идентификатора процесса на вашем устройстве. Следовательно, по умолчанию разные сеансы дают разные результаты моделирования.

Давайте снова попробуем сгенерировать 5 чисел из случайного нормального распределения, но в этот раз зафиксируем seed:

set.seed(123)
rnorm(5,0,1)

## [1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774

И еще разок:

set.seed(123)
rnorm(5,0,1)

## [1] -0.56047565 -0.23017749  1.55870831  0.07050839  0.12928774

Получилось! В качестве seed можно использовать любое число, главное, чтобы при репликации оно каждый раз было одинаковым.

Почитать про set.seed можно тут на русском и тут на английском.

0.5 4. Мини-симуляция

Симуляции – это “игрушечные” примеры, которые позволят нам на протяжении курса разбираться с разными методами. Игрушечные они потому, что за ними не стоят реальные данные. Данные для симуляций мы будем специальным образом заранее моделировать. Это бывает удобно, когда идеально подходящих данных нет, или они не лежат в открытом доступе. К тому же, живые данные часто могут быть зашумлены из-за других факторов, на которые нам не всегда будет удобно отвлекаться на занятиях. Но с живыми данными мы тоже обязательно будем работать! В том числе, они будут доступны вам в домашних заданиях :)

Теперь смоделируем небольшую симуляцию по мотивам изученого материала.

0.5.1 4.1. Подготовка данных

Давайте смоделируем гипотетическую ситуацию. Мы хотим оценить величину эффекта от использование сайта с расписанием cacs.ws на свободное время студента.

Предположим, что наша экспериментальная выборка состоит из 1000 человек:

N <- 1000

Для простоты у них будет всего две характеристики (\(X\) – возраст и \(Z\) – время в пути от дома до ЭФ), имеющих влияние на потенциальный исход (\(Y_0\) и \(Y_1\) – свободное время). При этом предположим, что реальный эффект воздействия равен \(\tau=15\) минутам. Будем считать, что реальная зависимость потенциального исхода от ковариатов и тритмента выглядит следующим образом:

• \(Y_0 = 240 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 0} + \varepsilon_0 = 120 - 3 \cdot X - Z + \varepsilon_0\) - потенциальный исход, если не было оказано воздействие
• \(Y_1 = 240 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 1} + \varepsilon_1 = 120 - 3 \cdot X - Z + 15 + \varepsilon_1\) - потенциальный исход, если было оказано воздействие

Сгенерируем \(X\) и \(Z\) случайным образом:

set.seed(123)
X <- runif(N, 18, 25)
Z <- rnorm(N, 60, 20)
# Посмотрим на наши данные, нарисуем сразу 2 графика -- диаграмму разброса и гистограмму
par(mfrow=c(2,2)) # читаем буквально как матрицы -- 2 строки и 2 столбца
plot(X)
hist(X)
plot(Z)
hist(Z)

Также создадим случайные ошибки для каждого из типов индивидов:

set.seed(123)
e0 <- rnorm(N/2, 0, 1)
e1 <- rnorm(N/2, 0, 1)
plot(e0, col="red")
points(e1, col='blue')

Вспомогательно пронумеруем наши наблюдения:

number <- 1:N

И соберем всё в общий датасет:

data <- data.frame(number = number, X=X, Z=Z) # сшиваем столбики в data frame
head(data, 10) # смотрим только первые несколько строк датасета

0.5.2 4.2. Рандомизация

Что касается рандомизации, то наша глобальная задача, используя разные функции придумать способ, как случайным образом разбить выборку на 2 группы. Существует огромное множество вариантов, как это можно сделать. Мы обсудим лишь несколько из них. На практике вы можете пользоваться как этими вариантами, так и придумать что-то своё.

0.5.2.1 4.2.1. Использование порядкового номера

Первый вариант – просто взять половину наблюдений по их порядку расположения в датасете. Пусть воздействие будет оказано только на первую половину нашей выборки.

T1 <- c(rep(1,N/2), rep(0,N/2))
mean(T1) # проверяем долю тритмента, должна получиться ровно половина

## [1] 0.5

data$T1 <- T1 # добавляем переменную в датасет

Тут мы воспользовались функцией c(...), которая позволяет объединить величины внутри неё в вектор, и функцией rep(x;times) которая повторяет величину \(x\) столько раз, сколько указано в величине \(times\).

Второй вариант, который практически дублирует первый, но немного отличается – мы можем взять каждое второе наблюдение:

T2 <- rep(0:1, N/2)
mean(T2) # проверяем долю тритмента, должна получиться ровно половина

## [1] 0.5

data$T2 <- T2 # добавляем переменную в датасет

0.5.2.2 4.2.2. Использование свойств распределений

Логично, что если мы хотим случайным образом присвоить наблюдениям разные группы, то мы можем начать с того, что присвоим им случайные числа, которые мы потом разными способами переведем в бинарный формат, который уже будет соотвествовать конкретной группе. Протестируем этот способ на примере трех распределений – нормального, равномерного и биномиального.

0.5.2.2.1 1) Равномерное распределение

set.seed(123)
V1 <- runif(N) # генерим вспомогательную переменную из равномерного распределения
T3 <- as.numeric(V1 < median(V1)) #генерим тритмент, отсекая половину выборки по медиане; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
mean(T3) # проверяем долю тритмента, должна получиться ровно половина

## [1] 0.5

data$T3 <- T3 # добавляем переменную в датасет

0.5.2.2.2 2) Нормальное распределение

set.seed(123)
V2 <- rnorm(N) # генерим вспомогательную переменную из нормального распределения
T4 <- as.numeric(V2>0) #генерим тритмент, отсекая половину выборки по знаку вспомогательной переменной; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
summary(T4) # проверяем долю тритмента, должна получиться примерно половина

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   0.505   1.000   1.000

data$T4 <- T4 # добавляем переменную в датасет

0.5.2.2.3 3) Биномиальное распределение

set.seed(123)
T5 <- rbinom(N, 1, 0.5) #генерим тритмент, отсекая половину выборки по знаку вспомогательной переменной; as.numeric используется, чтобы перейти от логического типа (true, false) к численному (1 и 0)
summary(T5) # проверяем долю тритмента, должна получиться примерно половина

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.493   1.000   1.000

data$T5 <- T5 # добавляем переменную в датасет

0.5.2.3 4.2.3. Использование готовых пакетов

0.5.2.3.1 1) Пакет experiment

library('experiment') 
set.seed(123)
rand <- randomize(data, group= c("Treat", "Control"))
T6 <- as.numeric(rand$treatment == "Treat") # преобразовываем переменную в численный формат 
data$T6 <- T6 # добавляем переменную в датасет
head(data, 10)

0.5.2.3.2 2) Пакет radomizr

library('randomizr') 
set.seed(123)
T7 <- complete_ra(N=N, m=N/2)
data$T7 <- T7 # добавляем переменную в датасет
head(data, 10)

0.5.2.4 4.2.4. Хэш-функции

Существует ряд ситуаций, когда подходы к рандомизации, которые мы разобрали выше, работают не очень хорошо. Например, в вашу выборку добавилось несколько наблюдений, которые по какой-то причине добавились не в конец вашего датасета. В этом случае рандомизация не будет воспроизводимой, а тритмент и контрольная группы будут различаться. Однако есть способ решить эту сложность.

Хэш-функция (hash function) – функция, осуществляющая преобразование массива входных данных произвольной длины в выходную битовую строку установленной длины, выполняемое определённым алгоритмом. Преобразование, производимое хэш-функцией, называется хэшированием.

Хэш-функции применяются в следующих случаях:

• при построении ассоциативных массивов;
• при поиске дубликатов в последовательностях наборов данных;
• при построении уникальных идентификаторов для наборов данных;
• при вычислении контрольных сумм от данных (сигнала) для последующего обнаружения в них ошибок (возникших случайно или внесённых намеренно), возникающих при хранении и/или передаче данных;
• при сохранении паролей в системах защиты в виде хэш-кода (для восстановления пароля - по хэш-коду требуется функция, являющаяся обратной по отношению к использованной хэш-функции);
• при выработке электронной подписи (на практике часто подписывается не само сообщение, а его “хэш-образ”);
• и др.

Для нас важно то, что с помощью хэш-функции мы сможем взаимно однозначно переводить данные формата строки в число.

Пример хэширования с помощью пакета digest:

library("digest")
hash <- digest('econometrics', algo="murmur32")
hash

## [1] "24b6de5a"

Пакет digest применяет хеш-функцию к произвольным объектам R. В пакете реализовано много разных алгоритмов преобразований, однако мы будем использовать алгоритм “murmur32” (подробнее почитать можно тут). Этот алгоритм совершенно не подходит для криптографических целей, но зато идеально подходит нам, поскольку он 32-битный, что позволяет нам перевести полученный хэш-код в целое число.

Прежде чем начать разбирать пример, немного отвлечемся на техническую полезную вещь – функцию sapply().

Про семейство apply функций рекомендую почитать подробнее на русском тут или тут, а на английском тут.

Прелесть функции sapply() состоит в том, что она позволяет нам избежать громоздких циклов при написании кода и ускорить вычисления благодаря “векторной ориентированности” языка R. Например, функция позволяет нам найти минимальное и максимальное значение для каждой из ковариат:

sapply(list(X,Z), min)

## [1] 18.003257  3.804506

sapply(list(X,Z), max)

## [1]  24.99583 127.80742

То есть функция проводит однотипную операцию (по сути вложенную функцию) над каждым элементом списка. Если кто-то хорошо владеет python, то аналогичной функцией там является map().

Вернемся к нашей симуляции. Изначально мы это все затеяли, чтобы сделать хорошую рандомизацию. Сначала хэшируем номера наших наблюдений. Так делать не очень хорошо, обычно, для хэширования данных используют ФИО и/или СНИЛС, но мы не будем ради этого громоздить генерацию еще одной переменной. Для иллюстрации просто воспользуемся номерами, это тоже сработает.

hashes <- sapply(data$number, function(x) {digest(x, algo="murmur32")}) # хэшируем строчки
hashes[1:20]

##  [1] "5e6216f3" "8b8b3789" "ad1bf356" "b2df92b2" "1d1b8e3a" "98d09cab"
##  [7] "2e366c5a" "4d110361" "91d6bfe0" "f568e543" "08bd74cd" "e979af4b"
## [13] "a82e5827" "485cbd08" "bc4e5bd8" "95f68430" "7180d744" "da414596"
## [19] "70aa2885" "afd0ffd2"

Далее нужно перевести получившиеся строчки в цифровой числовой. Функция strtoi() конвертирует строковое представление числа, которое хранится в строке, в длинное целое.

result <- strtoi(substring(hashes, 2), base=16) 
result[1:20]

##  [1] 241309427 193673097 219935574  48206514 219909690 147889323 238447706
##  [8] 219218785  30851040  90760515 146633933 158969675 137254951 140295432
## [15] 206461912 100041776  25220932 172049814  11151493 265355218

Далее, используя эти числа, нужно как-то разбить выборку на тритмент и контроль. Будем смотреть на последнюю цифру, если она меньше 5, то наблюдение попадет в тритмент, если больше – в контроль.

T8 <- result %% 10 < 5 # берем остаток от деления на 10 (получится последняя цифра) и сравниваем его с 5
data$T8 <- as.numeric(T8) # переводим из логического типа в числовой и добавляем в датасет
data$T8[1:20]

##  [1] 0 0 1 1 1 1 0 0 1 0 1 0 1 1 1 0 1 1 1 0

Проверим нашу рандомизацию на сбалансированность:

summary(data)

##      number             X               Z                 T1            T2     
##  Min.   :   1.0   Min.   :18.00   Min.   :  3.805   Min.   :0.0   Min.   :0.0  
##  1st Qu.: 250.8   1st Qu.:19.78   1st Qu.: 46.232   1st Qu.:0.0   1st Qu.:0.0  
##  Median : 500.5   Median :21.43   Median : 60.576   Median :0.5   Median :0.5  
##  Mean   : 500.5   Mean   :21.48   Mean   : 60.239   Mean   :0.5   Mean   :0.5  
##  3rd Qu.: 750.2   3rd Qu.:23.23   3rd Qu.: 73.101   3rd Qu.:1.0   3rd Qu.:1.0  
##  Max.   :1000.0   Max.   :25.00   Max.   :127.807   Max.   :1.0   Max.   :1.0  
##        T3            T4              T5              T6            T7     
##  Min.   :0.0   Min.   :0.000   Min.   :0.000   Min.   :0.0   Min.   :0.0  
##  1st Qu.:0.0   1st Qu.:0.000   1st Qu.:0.000   1st Qu.:0.0   1st Qu.:0.0  
##  Median :0.5   Median :1.000   Median :0.000   Median :0.5   Median :0.5  
##  Mean   :0.5   Mean   :0.505   Mean   :0.493   Mean   :0.5   Mean   :0.5  
##  3rd Qu.:1.0   3rd Qu.:1.000   3rd Qu.:1.000   3rd Qu.:1.0   3rd Qu.:1.0  
##  Max.   :1.0   Max.   :1.000   Max.   :1.000   Max.   :1.0   Max.   :1.0  
##        T8       
##  Min.   :0.000  
##  1st Qu.:0.000  
##  Median :1.000  
##  Mean   :0.508  
##  3rd Qu.:1.000  
##  Max.   :1.000

0.5.3 4.3. Реализация исходов

Посчитаем потенциальные исходы:

Y0 <- 240 - 3*X - Z + 15*0 + e0 # T=0
Y1 <- 240 - 3*X - Z + 15*1 + e1 # T=1
data$Y0 <- Y0 # добавляем переменные в датасет
data$Y1 <- Y1

А также наблюдаемые исходы:

data$Y <- T1*Y1 + (1-T1)*Y0

0.5.4 4.3. Считаем эффект

Когда данные подготовлены, мы забываем, что что-то о них знали :)

С данного момента мы действуем согласно предпосылке, что мы знаем только \(Y\), \(X\) и \(T\).

Поскольку мы знаем, что тритмент распределялся среди наблюдений случайно, мы можем сделать вывод, что предпосылка о независимости воздействия выполнена. Также мы знаем, что данные устроены так, что SUTVA тоже выполнена. Следовательно, можем использовать обычную разницу в средних, чтобы оценить эффнект воздействия:

\(\widehat{ATE} = \overline{Y_1} - \overline{Y_0}\)

ATE_hat <- mean(data$Y[T1==1]) - mean(data$Y[T1==0])
ATE_hat

## [1] 15.61754

Из курса ЭКМ-2 помним, что то же самое можно было бы получить, с помощью обычного парного МНК:

model1 <- lm(Y ~ T1, data=data)
summary(model1)

## 
## Call:
## lm(formula = Y ~ T1, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -76.984 -14.469  -0.314  13.235  60.998 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 115.0259     0.9094  126.49   <2e-16 ***
## T1           15.6175     1.2861   12.14   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 20.33 on 998 degrees of freedom
## Multiple R-squared:  0.1287, Adjusted R-squared:  0.1279 
## F-statistic: 147.5 on 1 and 998 DF,  p-value: < 2.2e-16

Если обе предпосылки выполнены, то оценка должна быть несмещенной даже без контрольных переменных (ковариат). Попробуем их добавить и сравним оценки эффекта:

model2 <- lm(Y ~ T1 + X + Z, data=data)
summary(model2)

## 
## Call:
## lm(formula = Y ~ T1 + X + Z, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9500 -0.4804 -0.0218  0.4296  3.2849 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 238.445186   0.296642   803.8   <2e-16 ***
## T1           14.978356   0.052761   283.9   <2e-16 ***
## X            -3.001684   0.013116  -228.9   <2e-16 ***
## Z            -0.973141   0.001318  -738.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8341 on 996 degrees of freedom
## Multiple R-squared:  0.9985, Adjusted R-squared:  0.9985 
## F-statistic: 2.266e+05 on 3 and 996 DF,  p-value: < 2.2e-16

Видим, что результаты устойчивы, величина эффекта практически не изменилась, но наша первая оценка все-таки была слегка завышеная, значит небольшое смещение было. Проверим контрольную и тритмент группу на “похожесть”, то есть проведем баланс ковариатов.

Далее в качестве упраженения посчитаем средние эффекты на разных группах:

ATT <- mean(data$Y1[data$T1 == 1] - data$Y0[data$T1 == 1])
ATT

## [1] 14.96307

ATnT <- mean(data$Y1[data$T1 == 0] - data$Y0[data$T1 == 0])
ATnT

## [1] 14.96307

Получили, то что и должны были получить – ATT = ATnT. Если мы вернемся в раздел, где мы подготовливали данные, то увидим, что в двух потенциальных исходах мы “зашифровали” одинаковый эффект:

• \(Y_0 = 120 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 0} + \varepsilon_0 = 120 - 3 \cdot X - Z + \varepsilon_0\)
• \(Y_1 = 120 - 3 \cdot X - Z + 15 \cdot \underbrace{T}_{= 1} + \varepsilon_1 = 120 - 3 \cdot X - Z + 15 + \varepsilon_1\)

Содержательно это иллюстрирует то, что эффект воздействия гомогенен, то есть \(\text{Heterogeneous treatment effect bias} = (1-\pi)(ATT - ATnT) = 0\).

0.5.5 4.4. Баланс ковариатов

Теперь проверим насколько наша рандомизация оказалась хорошей. Сравним средние значения ковариат в двух группах.

0.5.5.1 4.4.1. Сравнение средних

mean(data$X[T1==1])-mean(data$X[T1==0])

## [1] -0.02791759

mean(data$Z[T1==1])-mean(data$Z[T1==0])

## [1] -0.5707093

Это не очень информативный способ, но если мы вспомним, что порядок значений ковариат измерялся десятками, то можно сделать вывод, что группы в среднем похожи.

0.5.5.2 4.4.2. Тест Стьюдента

А теперь проведем тест на разницу в средних с помощью формального t-теста Стьюдента:

• \(H_0:\) среднее значение параметра в группах одинаковое, группы не различаются
• \(H_1:\) среднее значение параметра в группах разное, группы различаются

t.test(data$X[T1==1], data$X[T1==0])

## 
##  Welch Two Sample t-test
## 
## data:  data$X[T1 == 1] and data$X[T1 == 0]
## t = -0.21924, df = 997.51, p-value = 0.8265
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2777933  0.2219581
## sample estimates:
## mean of x mean of y 
##  21.46699  21.49490

t.test(data$Z[T1==1], data$Z[T1==0])

## 
##  Welch Two Sample t-test
## 
## data:  data$Z[T1 == 1] and data$Z[T1 == 0]
## t = -0.45032, df = 997.68, p-value = 0.6526
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.057656  1.916238
## sample estimates:
## mean of x mean of y 
##  59.95331  60.52402

Согласно значению p-value мы принимаем \(H_0\) и делаем вывод, что группы одинаковые.

0.5.5.3 4.4.2. Использование готовых пакетов

library("tableone")
table1 <- CreateTableOne(vars=c("X", "Z"), strata="T1", data=data, test=TRUE)
print(table1)

##                Stratified by T1
##                 0             1             p      test
##   n               500           500                    
##   X (mean (SD)) 21.49 (2.04)  21.47 (1.99)   0.827     
##   Z (mean (SD)) 60.52 (19.86) 59.95 (20.22)  0.653

0.5.6 4.5 Функции

Можно ли было бы сделать все то же самое за меньшее число строчек кода и время? Конечно, можно. Более того, если бы похожую процедуру вам пришлось бы по каким-то причинам проводить несколько раз, например, когда у вас несколько экспериментов, функция для вас просто была бы незаменимым помощником. Попробуем переписать все то же самое в несколько функций.

1. Функция для генерации данных

generate_data <- function(N){
  N=N
  X <- runif(N, 18, 25)
  Z <- rnorm(N, 60, 20)
  e0 <- rnorm(N/2, 0, 1)
  e1 <- rnorm(N/2, 0, 1)
  number <- 1:N
  Y0 <- 240 - 3*X - Z + 15*0 + e0 
  Y1 <- 240 - 3*X - Z + 15*1 + e1 
  data <- data.frame(number = number, X=X, Z=Z, Y0=Y0, Y1=Y1)
}

2. Функция для рандомизации

randomization <- function(data, type){
  N <- nrow(data)
  if (type=='in order'){
    T <- c(rep(1,N/2), rep(0,N/2))
    data$T <- T
  } else if (type=='by turns'){
    T <- rep(0:1, N/2)
    data$T <- T
  } else if (type=='unif'){
    V <- runif(N)
    T <- as.numeric(V1 < median(V1))
    data$T <- T
  } else if (type=='norm'){
    V <- rnorm(N) 
    T <- as.numeric(V2>0)
    data$T <- T
  } else if (type=='binom'){
    T <- rbinom(N, 1, 0.5)
  }
  return(data)
}

3. Функция для расчета наблюдаемого исхода

outcome <- function(data){
  observed_data <- data
  observed_data$Y <- data$T*observed_data$Y1 + (1 - data$T)*observed_data$Y0
  return(observed_data)
}

4. Функция для оценки эффекта

estimate <- function(data, response_variable) {
    m <- mean(data[data$T==1, response_variable])-mean(data[data$T==0, response_variable])
    return(m)
}

5. Функция для баланса ковариатов

balance <- function(data, covariates, type){
  if (type=='mean'){
    for (i in 1:length(covariates)){
      cov <- covariates[i]
      m <- mean(data[data$T==1, cov])-mean(data[data$T==0, cov])
      print(cov)
      print(m)
    }
  } else if (type=='t.test'){
    for (i in 1:length(covariates)){
      cov <- covariates[i]
      t <- t.test(data[data$T==1, cov], data[data$T==1, cov])
      print(cov)
      print(t)
    }
  }
}

А теперь посмотрим, сколько строчек займет наша симуляция с использованием функций:

set.seed(123)
data1 <- generate_data(1000)
set.seed(123)
data1 <- randomization(data=data1, type='in order')
data1 <- outcome(data=data1)
tail(data1,10)

estimate(data=data1, response_variable = 'Y')

## [1] 15.53169

balance(data=data1, covariates = c('X','Z'), type = 'mean')

## [1] "X"
## [1] -0.02791759
## [1] "Z"
## [1] -0.5707093

0.6 5. Тестирование гипотез

Ваш исследовательский вопрос может быть таким, что вам интересно оценить воздействия разных типов тритмента, то есть у вас есть несколько экспериментальных групп и одна контрольная. При такой постановке мы хотим проверить не одну, а сразу много статистических гипотез о различии в группах. При проверке любой гипотезы существует вероятность совершить ошибку первого рода (отклонить нулевую гипотезу, если она верна = обнаружить эффект, которого нет). Особенность множественного тестирования гипотез состоит в том, что чем больше гипотез мы проверяем на одних и тех же данных, тем больше будет вероятность допустить как минимум одну ошибку первого рода – эффект множественных сравнений (multiple comparisons/testing).

Рассмотрим это на примере. Предположим, что у нас есть 3 группы (A, B и С), в которых мы хотим сравнить среднее значение переменной интереса. Как и ранее, мы будем использовать t-тест Стьюдента. Если мы получили достаточно большое значение t-статистики такое, что p-value < 0.01, то мы отклоняем нулевую гипотезу и заключаем, что группы статистически различаются по переменной интереса. Отсечка p-value < 0.01 значит, что вероятность ошибочного вывода о различии между групповыми средними не превышает 0.01. Это будет работать именно так, когда у нас всего две группы, но в случае множественного тестирования вероятность будет больше 1%.

Выполняя тест Стьюдента, исследователь проверяет нулевую гипотезу об отсутствии разницы между двумя группами. Сравнивая группы A и В, он может ошибиться с вероятностью 1%, В и С – 1%, А и С – тоже 1%. Соответственно, вероятность ошибиться хотя бы в одном из этих трех сравнений составит:

\(P = 1 - \left(1-\alpha \right)^n = 1 - 0.99^3 \approx 0.03 > 0.01\)

Если бы групп было бы 5:

\(P = 1 - \left(1-\alpha \right)^n = 1 - 0.99^{10} \approx 0.1 > 0.01\)

К счастью, существует несколько методов, позволяющих преодолеть эту сложность.

0.6.1 5.1. Методы решения проблемы множественного тестирования гипотез

0.6.1.1 5.1.1. Групповая вероятность ошибки первого рода (family-wise error rate)

Предположим, что мы проверяем \(n\) гипотез. Для каждой гипотезы мы будем проводить тест Стьюдента. Результаты наших тестов можно обобщить следующим образом:

	Число принятых гипотез	Число отвергнутых гипотез	Всего
Число верных гипотез	\(U\) – Число безошибочно принятых гипотез (true negatives)	\(V\) – Число ошибочно отвергнутых гипотез (false positives) – ошибка первого рода	\(U+V=n^*\) – Число верных нулевых гипотез (true null hypotheses)
Число неверных гипотез	\(T\) – Число ошибочно принятых гипотез (false negatives) – ошибка второго рода)	\(S\) – Число безошибочно отвергнутых гипотез	\(T+S=n-n^*\) – Число истинных альтернативных гипотез (true alternative hypotheses)
Всего	\(U+T=W\) – Общее число принятых гипотез	\(V+S=R\) – Общее число отвергнутых гипотез

В теории всего существует \(n^*\) верных нулевых гипотез. В результате наших тестов мы ошибочно отвергаем \(V\) гипотез и верно принимаем остальные \(U\) гипотез. Также существует \(n−n^*\) альтернативных гипотез, из которых \(S\) гипотез безошибочно отвергаются, а \(T\) гипотез – ошибочно принимаются. Важно, что общие количества отвергнутых и принятых гипотез (\(R\) и \(W\)), а следовательно, и суммарное число гипотез \(n\) нам известны, тогда как остальные величины (\(n^*\), \(U\), \(V\), \(T\) и \(S\)) мы не наблюдаем.

При одновременной проверке семейства статистических гипотез мы хотим, чтобы количество наших ошибок (\(V\) и \(T\)) было минимальным. Традиционно исследователи пытаются минимизировать величину ошибочно отвергнутых гипотез \(V\). Это вполне логично, поскольку ложно отвергнутая нулевая гипотеза грозит нам ложноположительным найденным эффектом, которого реально может не быть.

Если \(V \geq 1\), мы совершаем как минимум одну ошибку первого рода. Вероятность допущения такой ошибки при множественной проверке гипотез называют групповой вероятностью ошибки (familywise error rate, FWER или experiment-wise error rate). По определению, \(FWER = P(V \geq 1)\). Соответственно, когда мы говорим, что хотим контролировать групповую вероятность ошибки на определенном уровне значимости \(\alpha\), мы подразумеваем, что должно выполняться неравенство \(FWER \leq \alpha\).

Ниже мы обсудим методы, которые позволяют это делать.

0.6.1.1.1 1) Коррекция Бонферрони

Вернемся к нашему примеру, когда мы сравнили 3 группы A, B и C с помощью t-теста. Предположим, что мы получили следующие Р-значения: 0.001, 0.01 и 0.04.

Как было сказано вышем, мы хотим, чтобы групповая вероятность ошибки была не больше уровня значимости \(FWER \leq \alpha\). Согласно методу Бонферрони, мы должны сравнить каждое из полученных p-значений не с \(\alpha\), а с \(\frac{\alpha}{n}\), где \(n\) – число проверяемых гипотез. Деление исходного уровня значимости \(\alpha\) на \(n\) – это и есть поправка Бонферрони. В рассматриваемом примере каждое из полученных p-значений необходимо было бы сравнить с \(\frac{\alpha}{n}\), например, с \(\frac{0.01}{3}\approx 0.0033\).

• \(p-value_1=0.001 < \alpha_{adjusted}=0.0033\) – гипотеза отклонена
• \(p-value_1=0.01 > \alpha_{adjusted}=0.0033\) – гипотеза принята
• \(p-value_1=0.04 > \alpha_{adjusted}=0.0033\) – гипотеза принята

Вместо деления уровня значимости на число гипотез, мы могли бы умножить каждое p-значение на это число и получить точно такие же выводы (эта эквивалентная процедура реалирована в R):

• \(p-value_{1,adjusted} = 0.001 \cdot 3 = 0.003 < \alpha = 0.01\) – гипотеза отклонена
• \(p-value_{1,adjusted} = 0.01 \cdot 3 = 0.03 > \alpha = 0.01\) – гипотеза принята
• \(p-value_{1,adjusted} = 0.04 \cdot 3 = 0.12 > \alpha = 0.01\) – гипотеза принята

Иногда при домножении p-значений результат может получиться больше единицы. Из теории вероятностей мы знаем, что вероятность не может быть больше одного, поэтому в таких случаях p-значение принимают равным за единицу.

Различные виды коррекций p-значений представлены в функции p.adjust(), выбрать тип коррекции можно с помощью аргумента method. В этой функции используется домножение исходных p-значений на количество тестируемых гипотез, а не корректировка уровня значимости.

Проверим наши рассчеты:

p.adjust(c(0.001, 0.01, 0.04), method = "bonferroni")

## [1] 0.003 0.030 0.120

Можно на выходе сразу получить выводы относительно гипотез при \(\alpha =5%\):

alpha <- 0.05
p.adjust(c(0.001, 0.01, 0.04), method = "bonferroni") < alpha # отклоняем H_0 (есть эффект)? 

## [1]  TRUE  TRUE FALSE

Важно помнить об уязвимости коррекции Бонферрони – с ростом числа гипотез мощность метода уменьшается. Чем больше гипотез мы хотим проверить, тем сложнее нам будет их отвергать (даже если они реально должны быть отвергнуты). Например, для 5 групп (10 гипотез), применение поправки Бонферрони привело бы к снижению исходного уровня значимости до 0.01/10 = 0.001. Соответственно, для отклонения гипотез, соответствующие p-значения должны быть меньше 0.001, а это довольно жесткая отсечка. Из этого делаем вывод, что при большом числе гипотез коррекцию Бонферрони лучше не использовать.

0.6.1.1.2 2) Низходящая процедура Хольма (Хольма-Бонферрони)

Метод Хольма позволяет побороть недостатки метода Бонферрони. Он устроен следующим образом:

• Сначала p-значения сортируются по возрастанию \(\displaystyle{p-value_1 \leq p-value_2 \leq ... \leq p-value_n}\).
• Затем проверяется условие для первого из p-значений: \(\displaystyle{p-value_1 \geq \frac{\alpha}{n-i+1}=\frac{\alpha}{n}}\), если условие выполнено, то все нулевые гипотезы принимаются, и процедура останавливается, иначе первая из гипотез отвергается, и начинается следующий шаг.
• На следующем шаге проверяется условие \(\displaystyle{p-value_2 \geq \frac{\alpha}{n-i+1}=\frac{\alpha}{n-1}}\), если условние выполнено, то все гипотезы, начиная со второй, принимаются, иначе первые две гипотезы отклоняются и начинается следующий шаг.
• На последнем шаге проверяется условие вида \(\displaystyle{p-value_n \geq \frac{\alpha}{n-n+1}}\), если оно выполнено, то последняя гипотеза принимается, если нет – отклоняется, на этом процедура заканчивается.

Метод Хольма называют нисходящей (step-down) процедурой. Он начинается с наименьшего p-значения в упорядоченном ряду и последовательно “спускается” вниз к более высоким значениям. На каждом шаге соответствующее значение \(p-value_i\) сравнивается со скорректированным уровнем значимости \(\displaystyle{\alpha_{adjusted}=\frac{\alpha}{n+i-1}}\). Аналогично коррекции Бонферрони можно вместо корректировки уровня значимости корректировать p-значения \(\displaystyle{p-value_{i,adjusted}=p-value_{i}\cdot(n-i+1)}\) (эта эквивалентная процедура реалирована в R). Возвращаясь к нашему примеру:

• \(p-value_{1,adjusted} = 0.001 \cdot (3-1+1) = 0.003 < \alpha = 0.01\) – гипотеза отклонена
• \(p-value_{1,adjusted} = 0.01 \cdot (3-2+1) = 0.02 > \alpha = 0.01\) – гипотеза принята
• \(p-value_{1,adjusted} = 0.04 \cdot (3-3+1) = 0.04 > \alpha = 0.01\) – гипотеза принята

А теперь проверим себя с помощью R:

p.adjust(c(0.001, 0.01, 0.04), method = "holm")

## [1] 0.003 0.020 0.040

И результаты проверки гипотез при \(\alpha =5%\):

alpha <- 0.05
p.adjust(c(0.001, 0.01, 0.04), method = "holm") < alpha # отклоняем H_0 (есть эффект)? 

## [1] TRUE TRUE TRUE

Про другие процедуры family-wise error rate можно почитать тут.

0.6.1.2 5.1.2. Средняя доля ложных отклонений

Рассмотренные выше FWER методы обеспечивают контроль над групповой вероятностью ошибки первого рода. Как мы выяснили, эти методы чересчур жестко работают, когда нужно одновременно проверить слишком много гипотез (падает статистическая мощность).Под “недостаточной мощностью” понимается сохранение многих нулевых гипотез, которые потенциально могут представлять исследовательский интерес и которые, соответственно, следовало бы отклонить. Недостаточная мощность традиционных процедур множественной проверки гипотез привела к разработке новых методов, например, метода Бенджамини-Хохберга.

0.6.1.2.1 1) Восходящая процедура Бенджамини — Хохберга

Для преодоления недостаточной мощности FWER методов был предложен новый подход к проблеме множественных проверок статистических гипотез. Суть подхода заключается в том, что вместо контроля над групповой вероятностью ошибки первого рода выполняется контроль над ожидаемой долей ложных отклонений (false discovery rate, FDR) среди всех отклоненных гипотез.

В терминах таблицы выше эта ожидаемая доля может быть записана следующим образом: \(\displaystyle{FDR=\left(\frac{V}{R}\right)}\).

В отличие от уровня значимости \(\alpha\), каких-либо общепринятых значений FDR не существует. Многие исследователи по аналогии контролируют FDR на уровне 5%. Интерпретация порогового значения FDR очень проста: например, если в ходе анализа данных отклонено 1000 гипотез, то при q=0.10 ожидаемая доля ложно отклоненных гипотез не превысит 100.

В статье (Benjamini, Hochberg, 1995) описание процедуры контроля над FDR выглядит так:

• Сначала p-значения сортируются по возрастанию \(\displaystyle{p-value_1 \leq p-value_2 \leq ... \leq p-value_n}\).
• Находят максимальное значение \(k\) среди всех индексов \(i=1,...,n\), для которого \(p-value_i \leq \frac{i}{m}q\) выполняется неравенство
• Отклоняют все гипотезы \(H_i\) с индексами \(i=1,...,k\)

Эквивалентная процедура, реалированая в R отличается тем, что вместо нахождения максимального индекса \(k\), исходные p-значения корректируются следующим образом: \(q_i=\frac{p_in}{i}\).

В качестве примера рассмотрим следующий ряд из 15 упорядоченных по возрастанию p-значений (из оригинальной статьи Benjamini and Hochberg 1995):

p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095,  0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BH")

##  [1] 0.00150000 0.00300000 0.00950000 0.03562500 0.06030000 0.06385714
##  [7] 0.06385714 0.06450000 0.07650000 0.48600000 0.58118182 0.71487500
## [13] 0.75323077 0.81321429 1.00000000

И результаты проверки гипотез при \(\alpha =5%\):

alpha <- 0.05
p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095,  0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BH") < alpha # отклоняем H_0 (есть эффект)? 

##  [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE

Интерпретация этих Р-значений с поправкой (в большинстве литературных источников их называют q-значениями) такова:

• Допустим, что мы хотим контролировать долю ложно отклоненных гипотез на уровне FDR = 0.05
• Все гипотезы, q-значения которых \(q-value \leq 0.05\), отклоняются
• Среди всех этих отклоненных гипотез доля отклоненных по ошибке не превышает 5%

Коррекция Р-значений по методу Беньямини-Хохберга работает особенно хорошо в ситуациях, когда необходимо принять общее решение по какому-либо вопросу при наличии информации (=проверенных гипотез) по многим параметрам.

Следует помнить, что описанный здесь метод контроля над ожидаемой долей ложных отклонений предполагает, что все тесты, при помощи которых получают p-значения, независимы. На практике в большинстве случаев это условие выполняться не будет.

0.6.1.2.2 2) Восходящая процедура Бенджамини-Йекутили

Для преодоления ограничения независимости тестов при проверке гипотез в работе (Benjamini and Yekutieli 2001) был предложен усовершенствованный метод, учитывающий наличие корреляции между проверяемыми гипотезами.

Процедура Бенджамини-Йекутили очень похожа на процедуру Бенджамини-Хохберга. Основное отличие заключается во введении поправочной константы \(\displaystyle{c_n=\sum \limits_{i=1}^{n}\frac{1}{i}}\), далее аналогично:

• Сначала p-значения сортируются по возрастанию \(\displaystyle{p-value_1 \leq p-value_2 \leq ... \leq p-value_n}\).
• Находят максимальное значение \(k\) среди всех индексов \(i=1,...,n\), для которого \(p-value_i \leq \frac{i}{m} \frac{q}{c_n}\) выполняется неравенство
• Отклоняют все гипотезы \(H_i\) с индексами \(i=1,...,k\)

В R реализуется эквивалентная процедура:

Эквивалентная процедура, реалированая в R отличается тем, что вместо нахождения максимального индекса \(k\), исходные p-значения корректируются следующим образом: \(\displaystyle{q_i=\frac{p_i\cdot n\cdot c_n}{i}}\).

p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095,  0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BY")

##  [1] 0.004977343 0.009954687 0.031523175 0.118211908 0.200089208 0.211892623
##  [7] 0.211892623 0.214025770 0.253844518 1.000000000 1.000000000 1.000000000
## [13] 1.000000000 1.000000000 1.000000000

И результаты проверки гипотез при \(\alpha =5%\):

alpha <- 0.05
p.adjust(c(0.0001, 0.0004, 0.0019, 0.0095,  0.0201, 0.0278, 0.0298, 0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000), method = "BY") < alpha # отклоняем H_0 (есть эффект)? 

##  [1]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE

0.6.1.3 5.1.3. Симуляция и сравнение результатов

Сравним как работают разные методы:

alpha <- 0.05
n <- 50
set.seed(123)
x <- rnorm(n, mean = c(rep(0, n/2), rep(3, n/2)))
pval <- round(2*pnorm(sort(-abs(x))), 3)

default_bool <- pval < alpha

bonferroni_pval <- p.adjust(pval, method = "bonferroni")
bonferroni_bool <- p.adjust(pval, method = "bonferroni") < alpha # отклоняем H_0 (есть эффект)? 

holm_pval <- p.adjust(pval, method = "holm")
holm_bool <- p.adjust(pval, method = "holm") < alpha # отклоняем H_0 (есть эффект)? 

bh_pval <- p.adjust(pval, method = "BH") 
bh_bool <- p.adjust(pval, method = "BH") < alpha # отклоняем H_0 (есть эффект)? 

by_pval <- p.adjust(pval, method = "BY")
by_bool <- p.adjust(pval, method = "BY") < alpha # отклоняем H_0 (есть эффект)? 

methods <- cbind(default_bool, bonferroni_bool, holm_bool, bh_bool, by_bool)
colnames(methods) <- c('Без коррекции', 'Бонферрони', 'Хольм', 'Бенджамини-Хохберг', 'Бенджамини-Йекутили')
rownames(methods) <- c(1:n)
methods

##    Без коррекции Бонферрони Хольм Бенджамини-Хохберг Бенджамини-Йекутили
## 1           TRUE       TRUE  TRUE               TRUE                TRUE
## 2           TRUE       TRUE  TRUE               TRUE                TRUE
## 3           TRUE       TRUE  TRUE               TRUE                TRUE
## 4           TRUE       TRUE  TRUE               TRUE                TRUE
## 5           TRUE       TRUE  TRUE               TRUE                TRUE
## 6           TRUE       TRUE  TRUE               TRUE                TRUE
## 7           TRUE       TRUE  TRUE               TRUE                TRUE
## 8           TRUE       TRUE  TRUE               TRUE                TRUE
## 9           TRUE       TRUE  TRUE               TRUE                TRUE
## 10          TRUE       TRUE  TRUE               TRUE                TRUE
## 11          TRUE      FALSE  TRUE               TRUE                TRUE
## 12          TRUE      FALSE FALSE               TRUE                TRUE
## 13          TRUE      FALSE FALSE               TRUE               FALSE
## 14          TRUE      FALSE FALSE               TRUE               FALSE
## 15          TRUE      FALSE FALSE               TRUE               FALSE
## 16          TRUE      FALSE FALSE               TRUE               FALSE
## 17          TRUE      FALSE FALSE               TRUE               FALSE
## 18          TRUE      FALSE FALSE               TRUE               FALSE
## 19          TRUE      FALSE FALSE               TRUE               FALSE
## 20          TRUE      FALSE FALSE               TRUE               FALSE
## 21          TRUE      FALSE FALSE              FALSE               FALSE
## 22          TRUE      FALSE FALSE              FALSE               FALSE
## 23         FALSE      FALSE FALSE              FALSE               FALSE
## 24         FALSE      FALSE FALSE              FALSE               FALSE
## 25         FALSE      FALSE FALSE              FALSE               FALSE
## 26         FALSE      FALSE FALSE              FALSE               FALSE
## 27         FALSE      FALSE FALSE              FALSE               FALSE
## 28         FALSE      FALSE FALSE              FALSE               FALSE
## 29         FALSE      FALSE FALSE              FALSE               FALSE
## 30         FALSE      FALSE FALSE              FALSE               FALSE
## 31         FALSE      FALSE FALSE              FALSE               FALSE
## 32         FALSE      FALSE FALSE              FALSE               FALSE
## 33         FALSE      FALSE FALSE              FALSE               FALSE
## 34         FALSE      FALSE FALSE              FALSE               FALSE
## 35         FALSE      FALSE FALSE              FALSE               FALSE
## 36         FALSE      FALSE FALSE              FALSE               FALSE
## 37         FALSE      FALSE FALSE              FALSE               FALSE
## 38         FALSE      FALSE FALSE              FALSE               FALSE
## 39         FALSE      FALSE FALSE              FALSE               FALSE
## 40         FALSE      FALSE FALSE              FALSE               FALSE
## 41         FALSE      FALSE FALSE              FALSE               FALSE
## 42         FALSE      FALSE FALSE              FALSE               FALSE
## 43         FALSE      FALSE FALSE              FALSE               FALSE
## 44         FALSE      FALSE FALSE              FALSE               FALSE
## 45         FALSE      FALSE FALSE              FALSE               FALSE
## 46         FALSE      FALSE FALSE              FALSE               FALSE
## 47         FALSE      FALSE FALSE              FALSE               FALSE
## 48         FALSE      FALSE FALSE              FALSE               FALSE
## 49         FALSE      FALSE FALSE              FALSE               FALSE
## 50         FALSE      FALSE FALSE              FALSE               FALSE

plot(pval, bonferroni_pval, col = "orange", type="p", pch=1)
lines(pval, holm_pval, col="green", type="p", pch=1)
lines(pval, bh_pval, col="blue", type="p", pch=1)
lines(pval, by_pval, col="violet", type="p", pch=1)
abline(h=alpha, col="red")
abline(v=alpha, col="red")
legend(x=0.6, y=0.5, 
       legend=c('Бонферрони', 'Хольм', 'Бенджамини-Хохберг', 'Бенджамини-Йекутили', 'Уровень значимости'), 
       col=c("orange", "green", "blue", "violet", "red"), 
       bty = "n", pch=1)

0.7 6. Планирование эксперимента (НЕ ГОТОВО)

Ранее мы говорили, что при достаточном объеме выборки и хорошей рандомизации в среднем тритмент группа и контрольная группа окажутся похожимы друг на друга за исключением оказываемого на них воздействия, в этом случае в качестве хорошей оценки эффекта можно использовать простую разницу в средних исходах.

По разным причинам может получиться так, что тритмент группа и контрольная группа не похожи друг на друга.

Существует 2 способа как бороться с несбалансированностью групп:

• ex ante: “быть умным заранее”. Изначально проводим эксперимент так, чтобы несбалансированности не было и используем простую оценку
  • Престратификация
• ex post: берем такие данные, какие уже есть, но используем какие-то более умные и сложные оценки, сложнее обычной разницы в средних

0.7.1 6.1 Престратификация

0.7.2 6.2 Постстратификация

0.8 7. Контрольные переменные

0.9 8. Плохой контроль

Если у вас возникли какие-то вопросы, или вы нашли неточности в тексте, напишите мне об этом на почту annastavnychuk@gmail.com